---
title: "Analysis with reverse coding"
author: "Alisa Bedrov"
date: "3/6/2022"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(haven)
library(psych)
library(ggpubr)
library(apaTables)
library(tidyverse)
library(rstatix)
library(MplusAutomation)
library(here)
library(gt)
library(sjPlot)
library(DiagrammeR)
library(semPlot)
library(glue)
library(corrplot)
library(sjlabelled)

data <- read_csv("study3_dat.csv")
data$lfe <- 8 - data$lfe
data$cnsq <- 8-data$cnsq
data$aut <- 8-data$aut
data$intr <- 8-data$intr
data$dstn <- 8-data$dstn
```
### Data Screening
```{r}
data %>% select(-ID) %>% summary()
data %>% select(-ID) %>% describe()

data %>% shapiro_test(rum)
ggdensity(data$rum, fill="lightgray")
```

### EFA Analysis for Study 3 Data
```{r}
efa <- data %>% select(-ID, -burden)
efa_1 <- mplusObject(
  
  TITLE = "Model 1 - Default Rotation;",
  
  VARIABLE = "usevar = rum strs ef dif dstr adj lie avd dstn intr aut obl exp glt uncf rep lfe cnsq;",
  
  ANALYSIS = "type = efa 1 5;
  estimator = MLR;
  parallel=50;",
  
  MODEL = "",
  
  PLOT = "type = plot3;",
  
  OUTPUT = "sampstat;",
  
  usevariables = colnames(efa),
  rdata=efa
)

efa_1_fit <- mplusModeler(efa_1, 
                          dataout=("efa_1.dat"),
                          modelout=("efa_1.inp"),
                          check=TRUE, run = TRUE, hashfilename=FALSE)

model_fit <- LatexSummaryTable(efa_1_fit,
                               keepCols=c("Title", "Parameters", "LL", "ChiSqM_Value", "ChiSqM_DF", "ChiSqM_PValue", "RMSEA_Estimate", "RMSEA_90CI_LB", "RMSEA_90CI_UB", "CFI", "TLI", "SRMR")) %>%
  mutate(Title = c("1-Factor", "2-Factor", "3-Factor", "4-Factor", "5-Factor")) %>%
  mutate_at(vars(contains("RMSEA")), ~format(., nsmall = 3)) %>% unite(CI, RMSEA_90CI_LB:RMSEA_90CI_UB, sep=", ", remove=TRUE) %>% 
  mutate(CI=paste0("(", CI, ")")) %>% 
  unite (RMSEA, RMSEA_Estimate:CI, sep=" ", remove=TRUE)

model_fit %>% gt() %>% tab_header(title = md("**Table 1**"), subtitle=md("*Summary of Model Fit Indices*")) %>% cols_label(Title = "Model", Parameters = md("Par"), LL = md("*LL*"), ChiSqM_Value = md("CHi^2"), ChiSqM_PValue = md("*p-value*"), ChiSqM_DF = md("*df*"), RMSEA = "RMSEA (90% CI)") %>% tab_options(column_labels.font.weight="bold") %>% fmt(c(6), fns=function(x) ifelse(x<.001, "<.001", scales::number(x,accuracy=0.01)))

x <- list(EFA=efa_1_fit$results$gh5$efa$eigenvalues,
          Parallel=efa_1_fit$results$gh5$efa$parallel_average) %>% as_tibble()

plot_data <- cbind(Factor = paste0(1:nrow(x)), x) %>% 
  mutate(Factor = fct_inorder(Factor)) %>% 
  pivot_longer(EFA:Parallel,
               names_to = "Analysis",
               values_to = "Eigenvalues")

plot_data %>% ggplot(aes(y=Eigenvalues,
                         x=Factor,
                         group=Analysis,
                         color=Analysis)) +
  geom_point() + geom_line() + theme_minimal() +
  labs(title = "Figure 1: Parallel Eigenvalue Plot")

ggsave(("eigenvalue_elbow_rplot2.png"), dpi=300, height=5, width=7, units="in")


# Eigenvalues/Factor loadings
loadings_stdyx <- efa_1_fit$results$parameters$efa$f4$loadings$estimates %>% as_tibble() %>% rownames_to_column("Names") %>%   mutate(Names = str_to_lower(Names))

loadings_stdyx %>% gt() %>% 
  tab_header(title = md("**Table 2**"), 
             subtitle = md("*Summary of Factor Loadings: 4-Factor EFA Model*")) %>%
  cols_label('1'="F1", '2'="F2", '3'="F3", '4'="F4") %>% 
  tab_row_group(label="Factor 1: Daily Concealment", rows=1:8) %>% tab_row_group(label="Factor 2: Harm to Relationships", rows=9:11) %>% tab_row_group(label="Factor 3: Pull to Reveal", rows = 12:14) %>% tab_row_group(label="Factor 4: Anticipated Consequences for Self", rows = 15:18) %>%  row_group_order(groups=c("Factor 1: Daily Concealment", "Factor 2: Harm to Relationships", "Factor 3: Pull to Reveal", "Factor 4: Anticipated Consequences for Self")) %>% tab_style(style=list(cell_text(weight="bold")),locations = cells_body(columns="1", rows=1:8)) %>% tab_style(style=list(cell_text(weight="bold")),locations = cells_body(columns="2", rows=9:11)) %>% tab_style(style=list(cell_text(weight="bold")),locations = cells_body(columns="3", rows=12:14)) %>% tab_options(column_labels.font.weight="bold") %>% tab_style(style=list(cell_text(weight="bold")),locations = cells_body(columns="4", rows=15:18))
```

# APA style tables
### APA style tables
```{r}
apa.cor.table(efa, filename="figures2","Table1.doc", table.number=1)

apa <- function(x, title = " ") {
  gt(x) %>% 
    tab_options(table.border.top.color = "white",
                heading.title.font.size = px(16),
                column_labels.border.top.color = "black",
                column_labels.border.bottom.width = 3,
                column_labels.border.bottom.color = "black",
                table_body.border.bottom.color = "black",
                table.border.bottom.color = "white",
                table.width=pct(70),
                table.background.color = "white") %>% cols_align(align="center") %>% tab_style(style= list(
                  cell_borders(
                    sides=c("top","bottom"),
                    color="white",
                    weight=px(1)),
                  cell_text(align="center"),
                  cell_fill(color="white", alpha= NULL)),
                  locations=cells_body(columns=everything(),rows=everything())) %>% tab_header(title = html("<i>", title, "</i>")) %>% opt_table_font(font = "Times New Roman") %>% opt_align_table_header(align = "left")
}

fa_table <- loadings_stdyx %>% apa() %>% tab_header(title = md("**Table 2**"), subtitle = md("*Summary of Factor Loadings: 4-Factor EFA Model*")) %>% cols_label('1'="F1", '2'="F2", '3'="F3", '4'="F4") %>% tab_row_group(label="Factor 1: Daily Concealment", rows=c(1:8)) %>% tab_row_group(label="Factor 2: Harm to Relationships", rows=c(9:11)) %>% tab_row_group(label="Factor 3: Pull to Reveal", rows = c(12:14)) %>% tab_row_group(label = "Factor 4: Anticipated Consequences for Self", rows = c(15:18)) %>%  row_group_order(groups=c("Factor 1: Daily Concealment", "Factor 2: Harm to Relationships", "Factor 3: Pull to Reveal", "Factor 4: Anticipated Consequences for Self")) %>% tab_style(style=list(cell_text(weight="bold")),locations = cells_body(columns="1", rows=c(1:8))) %>% tab_style(style=list(cell_text(weight="bold")),locations = cells_body(columns="2", rows=c(9:11))) %>% tab_style(style=list(cell_text(weight="bold")),locations = cells_body(columns="3", rows=c(12:14))) %>% tab_style(style=list(cell_text(weight="bold")), locations=cells_body(columns="4", rows=c(15:18))) %>%  tab_options(column_labels.font.weight="bold")

fa_table

webshot::install_phantomjs(force=TRUE)

fa_table %>% gtsave("table2a.png")

fa_fit <- model_fit %>% apa() %>% tab_header(title = md("**Table 1**"), subtitle=md("*Summary of Model Fit Indices*")) %>% cols_label(Title = "Model", Parameters = md("Par"), LL = md("*LL*"), ChiSqM_Value = md("CHi^2"), ChiSqM_PValue = md("*p-value*"), ChiSqM_DF = md("*df*"), RMSEA = "RMSEA (90% CI)") %>% tab_options(column_labels.font.weight="bold") %>% fmt(c(6), fns=function(x) ifelse(x<.001, "<.001", scales::number(x,accuracy=0.01)))

fa_fit

fa_fit %>% gtsave("table3a.png")

```


### CFA Analysis with Study 4
```{r}
data2 <- read_csv("study4_dat.csv")
data2 %>% select(-ID) %>% describe()

data2$lfe <- 8 - data2$lfe
data2$cnsq <- 8 - data2$cnsq
data2$aut <- 8-data2$aut
data2$intr <- 8-data2$intr
data2$dstn <- 8-data2$dstn

cfa <- data2 %>% select(-ID, -burden)

grViz("digraph cfa_model {
      graph[layout=dot, overlap=true]
      
      node[shape=box]
      rum strs ef dif dstr adj lie avd dstn intr aut obl exp glt uncf rep lfe cnsq;
      
      node[shape=circle]
      factor1 factor2 factor3 factor4;
      
      edge[]
      factor1 ->{rum strs ef dif dstr adj lie avd}
      factor2 -> {dstn intr aut}
      factor3 -> {obl exp glt}
      factor4 -> {uncf rep lfe cnsq}
      
      edge[minlen=0]
     factor1 -> factor2 [arrowhead = none, tailport = 'n', headport = 'n']
      factor3 -> factor1 [arrowhead = none, tailport = 'n', headport = 'n']
      factor1 -> factor4 [arrowhead = none, tailport = 'n', headport = 'n']
      factor2 -> factor3 [arrowhead = none, tailport = 'n', headport = 'n']
      factor3 -> factor2 [arrowhead = none, tailport = 'n', headport = 'n']
      #factor2 -> factor4 [arrowhead = none, tailport = 'n', headport = 'n']
      #factor4 -> factor3 [arrowhead = none, tailport = 'n', headport = 'n']
      }")

cfa_m1 <- mplusObject(
  
  TITLE = "CFA - ULI",
  
  VARIABLE = "usevar = rum-cnsq;",
  
  ANALYSIS = "estimator = mlr;",
  
  MODEL = "factor1 by rum strs ef dif dstr adj lie avd;
  factor2 by dstn intr aut;
  factor3 by obl exp glt;
  factor4 by uncf rep lfe cnsq;",
  
  PLOT = "type = plot3;",
  
  OUTPUT = "sampstat standardized residual modindices (3.84);",
  
  usevariables = colnames(cfa),
  rdata = cfa)

cfa_m1_fit <- mplusModeler(cfa_m1,
                           dataout=("cfa_data.dat"),
                           modelout=("cfa_m1_uli.inp"),
                           check=TRUE, run = TRUE, hashfilename=FALSE)


out_uli <- readModels("cfa_m1_uli.out")
semPaths(out_uli, intercepts=FALSE)

model_fit_uli <- LatexSummaryTable(out_uli,
       keepCols=c("Title", "Parameters","LL", "ChiSqM_Value", "ChiSqM_DF", "ChiSqM_PValue", "RMSEA_Estimate", "RMSEA_90CI_LB", "RMSEA_90CI_UB", "CFI", "TLI", "SRMR")) %>% 
  mutate_at(vars(contains("RMSEA")), ~format(., nsmall=3)) %>% unite(CI, RMSEA_90CI_LB:RMSEA_90CI_UB, sep = ", ", remove=TRUE) %>% mutate(CI = paste0("(", CI, ")")) %>% unite(RMSEA, RMSEA_Estimate:CI, sep=" ", remove=TRUE)

model_fit_uli %>% gt() %>% tab_header(
  title = md("**Table 1**"),
  subtitle=md("*SUmmary of Model Fit Indices*")) %>% 
  cols_label(Title = "Model", Parameters = md("Par"),
             LL = md("*LL*"), ChiSqM_Value = md("CHi^2"),
             ChiSqM_PValue = md("*p-value*"), 
             ChiSqM_DF = md("*df*"),
             RMSEA = "RMSEA (90% CI)") %>% 
  tab_options(column_labels.font.weight = "bold") %>% 
  fmt(c(6), fns=function(x) ifelse(x<.001, "<.001", scales::number(x,accuracy=.01)))

var <- out_uli$parameters$stdyx.standardized %>% 
  filter(grepl(".BY", paramHeader)) %>% 
  select(param, est, se)

fac <- out_uli$parameters$stdyx.standardized %>% 
  filter(grepl(".WITH",paramHeader)) %>% 
  select(param, paramHeader, est, se) %>% 
  unite(param, paramHeader, param) %>% 
  mutate(param=str_replace(param, "FACTOR2.WITH_FACTOR1", "F1 with F2"),
         param=str_replace(param, "FACTOR3.WITH_FACTOR1", "F1 with F3"),
         param=str_replace(param, "FACTOR3.WITH_FACTOR2", "F2 with F3"),
         param=str_replace(param, "FACTOR4.WITH_FACTOR1", "F1 with F4"),
         param=str_replace(param, "FACTOR4.WITH_FACTOR3", "F3 with F4"),
         param=str_replace(param, "FACTOR4.WITH_FACTOR2", "F2 with F4"))

loadings_stdyx <- full_join(var,fac)

loadings_stdyx %>% gt() %>% tab_header(
  title = md("**Table 2**"),
  subtitle=md("*Standardized Factor Loadings and Factor Correlation Estimates*")) %>% 
  cols_label(param = "Item", est = "Loading", se = "SE") %>% 
  tab_row_group(label = "Factor 1", rows = 1:8) %>% 
  tab_row_group(label = "Factor 2", rows = 9:11) %>% 
  tab_row_group(label = "Factor 3", rows = 12:14) %>% 
  tab_row_group(label = "Factor 4", rows = 15:18) %>% 
  tab_row_group(label = "Factor Correlation", rows = 19:24) %>% 
  row_group_order(groups = c("Factor 1", "Factor 2", "Factor 3", "Factor 4", "Factor Correlation")) %>% 
  tab_options(column_labels.font.weight = "bold")

#Modification Indices
mod_ind <- out_uli$mod_indices %>% arrange(desc(MI))

mod_ind %>% gt() %>% tab_header(title=md("**Table 3**"),
     subtitle = md("*Model Modification Indices*")) %>% 
  cols_label(modV1 = "Variable 1", operator = "Operator",
             modV2 = "Variable 2", EPC = "E.P.C.", 
             Std_EPC = "Std E.P.C.", 
             StdYX_EPC = "StdYX E.P.C.") %>% 
  tab_options(column_labels.font.weight = "bold")


cfa_m3 <- mplusObject(
  
  TITLE = "CFA - ULI MOD",
  
  VARIABLE = "usevar = rum-cnsq;",
  
  ANALYSIS = "estimator = mlr;",
  
  MODEL = "factor1 by rum strs ef dif dstr adj lie avd;
  factor2 by dstn intr aut;
  factor3 by obl exp glt;
  factor4 by uncf rep lfe cnsq;
  cnsq with lfe",
  
  PLOT = "type = plot3;",
  
  OUTPUT = "sampstat standardized residual modindices (3.84);",
  
  usevariables = colnames(cfa),
  rdata = cfa)

cfa_m3_fit <- mplusModeler(cfa_m3,
                           dataout=("cfa_data.dat"),
                           modelout=("cfa_m3_uli.inp"),
                           check=TRUE, run = TRUE, hashfilename=FALSE)



out_mod <- readModels("cfa_m3_uli.out")


model_fit_mod <- LatexSummaryTable(out_mod,
       keepCols=c("Title", "Parameters","LL", "ChiSqM_Value", "ChiSqM_DF", "ChiSqM_PValue", "RMSEA_Estimate", "RMSEA_90CI_LB", "RMSEA_90CI_UB", "CFI", "TLI", "SRMR")) %>% 
  mutate_at(vars(contains("RMSEA")), ~format(., nsmall=3)) %>% unite(CI, RMSEA_90CI_LB:RMSEA_90CI_UB, sep = ", ", remove=TRUE) %>% mutate(CI = paste0("(", CI, ")")) %>% unite(RMSEA, RMSEA_Estimate:CI, sep=" ", remove=TRUE)

model_fit_mod %>% gt() %>% tab_header(
  title = md("**Table 1**"),
  subtitle=md("*Summary of Model Fit Indices*")) %>% 
  cols_label(Title = "Model", Parameters = md("Par"),
             LL = md("*LL*"), ChiSqM_Value = md("CHi^2"),
             ChiSqM_PValue = md("*p-value*"), 
             ChiSqM_DF = md("*df*"),
             RMSEA = "RMSEA (90% CI)") %>% 
  tab_options(column_labels.font.weight = "bold") %>% 
  fmt(c(6), fns=function(x) ifelse(x<.001, "<.001", scales::number(x,accuracy=.01)))

allFit <- full_join(model_fit_uli, model_fit_mod)

fa_fit2 <- allFit %>% apa() %>% tab_header(title = md("**Table 2**"),
    subtitle = md("*Summary of Model Fit Indices*"))  %>% 
  cols_label(Title = "Model", Parameters = md("Par"),
             LL = md("*LL*"), ChiSqM_Value = md("CHi^2"),
             ChiSqM_PValue = md("*p-value*"), 
             ChiSqM_DF = md("*df*"),
             RMSEA = "RMSEA (90% CI)") %>% 
  tab_options(column_labels.font.weight = "bold") %>% 
  fmt(c(6), fns=function(x) ifelse(x<.001, "<.001", scales::number(x,accuracy=.01)))

fa_fit2

fa_fit2 %>% gtsave("table4a.png")

# Mod INdices Part 2
mod_ind2 <- out_mod$mod_indices %>% arrange(desc(MI))

mod_ind2 %>% gt() %>% tab_header(title=md("**Table 3**"),
     subtitle = md("*Model Modification Indices*")) %>% 
  cols_label(modV1 = "Variable 1", operator = "Operator",
             modV2 = "Variable 2", EPC = "E.P.C.", 
             Std_EPC = "Std E.P.C.", 
             StdYX_EPC = "StdYX E.P.C.") %>% 
  tab_options(column_labels.font.weight = "bold")

cfa_m4 <- mplusObject(
  
  TITLE = "CFA - ULI MOD2",
  
  VARIABLE = "usevar = rum-cnsq;",
  
  ANALYSIS = "estimator = mlr;",
  
  MODEL = "factor1 by rum strs ef dif dstr adj lie avd;
  factor2 by dstn intr aut;
  factor3 by obl exp glt;
  factor4 by uncf rep lfe cnsq;
  cnsq with lfe;
  avd with adj",
  
  PLOT = "type = plot3;",
  
  OUTPUT = "sampstat standardized residual modindices (3.84);",
  
  usevariables = colnames(cfa),
  rdata = cfa)

cfa_m4_fit <- mplusModeler(cfa_m4,
                           dataout=("cfa_data.dat"),
                           modelout=("cfa_m4_uli.inp"),
                           check=TRUE, run = TRUE, hashfilename=FALSE)


out_mod2 <- readModels("cfa_m4_uli.out")

var2 <- out_mod2$parameters$stdyx.standardized %>% 
  filter(grepl(".BY", paramHeader)) %>% 
  select(param, est, se)

fac2 <- out_mod2$parameters$stdyx.standardized %>% 
  filter(grepl(".WITH",paramHeader)) %>% 
  select(param, paramHeader, est, se) %>% 
  unite(param, paramHeader, param) %>% 
  mutate(param=str_replace(param, "FACTOR2.WITH_FACTOR1", "F1 with F2"),
         param=str_replace(param, "FACTOR3.WITH_FACTOR1", "F1 with F3"),
         param=str_replace(param, "FACTOR3.WITH_FACTOR2", "F2 with F3"),
         param=str_replace(param, "FACTOR4.WITH_FACTOR1", "F1 with F4"),
         param=str_replace(param, "FACTOR4.WITH_FACTOR3", "F3 with F4"),
         param=str_replace(param, "FACTOR4.WITH_FACTOR2", "F2 with F4"))

loadings_stdyx2 <- full_join(var2,fac2)

loadings_stdyx2 %>% gt() %>% tab_header(
  title = md("**Table 2**"),
  subtitle=md("*Standardized Factor Loadings and Factor Correlation Estimates*")) %>% 
  cols_label(param = "Item", est = "Loading", se = "SE") %>% 
  tab_row_group(label = "Factor 1", rows = 1:8) %>% 
  tab_row_group(label = "Factor 2", rows = 9:11) %>% 
  tab_row_group(label = "Factor 3", rows = 12:14) %>% 
  tab_row_group(label = "Factor 4", rows = 15:18) %>% 
  tab_row_group(label = "Factor Correlation", rows = 19:24) %>% 
  tab_row_group(label = "Residual Correlations", rows=25:26) %>% 
  row_group_order(groups = c("Factor 1", "Factor 2", "Factor 3", "Factor 4", "Factor Correlation", "Residual Correlations")) %>% 
  tab_options(column_labels.font.weight = "bold")


model_fit_mod2 <- LatexSummaryTable(out_mod2,
       keepCols=c("Title", "Parameters","LL", "ChiSqM_Value", "ChiSqM_DF", "ChiSqM_PValue", "RMSEA_Estimate", "RMSEA_90CI_LB", "RMSEA_90CI_UB", "CFI", "TLI", "SRMR")) %>% 
  mutate_at(vars(contains("RMSEA")), ~format(., nsmall=3)) %>% unite(CI, RMSEA_90CI_LB:RMSEA_90CI_UB, sep = ", ", remove=TRUE) %>% mutate(CI = paste0("(", CI, ")")) %>% unite(RMSEA, RMSEA_Estimate:CI, sep=" ", remove=TRUE)

model_fit_mod2 %>% gt() %>% tab_header(
  title = md("**Table 1**"),
  subtitle=md("*Summary of Model Fit Indices*")) %>% 
  cols_label(Title = "Model", Parameters = md("Par"),
             LL = md("*LL*"), ChiSqM_Value = md("CHi^2"),
             ChiSqM_PValue = md("*p-value*"), 
             ChiSqM_DF = md("*df*"),
             RMSEA = "RMSEA (90% CI)") %>% 
  tab_options(column_labels.font.weight = "bold") %>% 
  fmt(c(6), fns=function(x) ifelse(x<.001, "<.001", scales::number(x,accuracy=.01)))

allFit2 <- full_join(allFit, model_fit_mod2)

fa_fit3 <- allFit2 %>% apa() %>% tab_header(title = md("**Table 3**"),
    subtitle = md("*Summary of Model Fit Indices*"))  %>% 
  cols_label(Title = "Model", Parameters = md("Par"),
             LL = md("*LL*"), ChiSqM_Value = md("CHi^2"),
             ChiSqM_PValue = md("*p-value*"), 
             ChiSqM_DF = md("*df*"),
             RMSEA = "RMSEA (90% CI)") %>% 
  tab_options(column_labels.font.weight = "bold") %>% 
  fmt(c(6), fns=function(x) ifelse(x<.001, "<.001", scales::number(x,accuracy=.01)))

fa_fit3

fa_fit3 %>% gtsave("table5.png")

mod_ind3 <- out_mod2$mod_indices %>% arrange(desc(MI))

mod_ind3 %>% gt() %>% tab_header(title=md("**Table 3**"),
     subtitle = md("*Model Modification Indices*")) %>% 
  cols_label(modV1 = "Variable 1", operator = "Operator",
             modV2 = "Variable 2", EPC = "E.P.C.", 
             Std_EPC = "Std E.P.C.", 
             StdYX_EPC = "StdYX E.P.C.") %>% 
  tab_options(column_labels.font.weight = "bold")

```

### Chi Square COmparisons
```{r}
c0 <- out_uli$summaries$ChiSqM_ScalingCorrection
c1 <- out_mod2$summaries$ChiSqM_ScalingCorrection

d0 <- out_uli$summaries$ChiSqM_DF
d1 <- out_mod2$summaries$ChiSqM_DF

SB0 <- out_uli$summaries$ChiSqM_Value
SB1 <- out_mod2$summaries$ChiSqM_Value

cd <- (((d0*c0)-(d1*c1))/(d0-d1))

TRd <- (((SB0*c0)-(SB1*c1))/(cd))
TRd2 <- abs(TRd)

df <- abs(d0-d1)

p_diff <- stats::pchisq(q=TRd2, df=df, lower.tail=FALSE)

p_diff


c0 <- out_uli$summaries$ChiSqM_ScalingCorrection
c2 <- out_mod$summaries$ChiSqM_ScalingCorrection

d0 <- out_uli$summaries$ChiSqM_DF
d2 <- out_mod$summaries$ChiSqM_DF

SB0 <- out_uli$summaries$ChiSqM_Value
SB2 <- out_mod$summaries$ChiSqM_Value

cd2 <- (((d0*c0)-(d2*c2))/(d0-d2))

TRda <- (((SB0*c0)-(SB2*c2))/(cd2))
TRd2a <- abs(TRda)

df2 <- abs(d0-d2)

p_diff2 <- stats::pchisq(q=TRd2a, df=df2, lower.tail=FALSE)

p_diff2


cd3 <- (((d1*c1)-(d2*c2))/(d1-d2))
TRD3 <- (((SB1*c1)-(SB2*c2))/(cd3))
TRD3a <- abs(TRD3)
df3 <- abs(d1-d2)
p_diff3 <- pchisq(TRD3a, df3, lower.tail=FALSE)
p_diff3

```

### Strictly Positive Satorra-Bentler Tests
```{r}
#cfa_m1 is original model, cfa_m3 has 1 mod ind, cfa_m4 has 2 mod ind
#re-run cfa_m1 model with svalues in output
cfa_m1SB <- mplusObject(
  
  TITLE = "CFA - ULI for SB Pos",
  
  VARIABLE = "usevar = rum-cnsq;",
  
  ANALYSIS = "estimator = mlr;",
  
  MODEL = "factor1 by rum strs ef dif dstr adj lie avd;
  factor2 by dstn intr aut;
  factor3 by obl exp glt;
  factor4 by uncf rep lfe cnsq;",
  
  PLOT = "type = plot3;",
  
  OUTPUT = "svalues",
  
  usevariables = colnames(cfa),
  rdata = cfa)

cfa_m1SB_fit <- mplusModeler(cfa_m1SB,
                           dataout=("cfaSB_data.dat"),
                           modelout=("cfa_m1SB_uli.inp"),
                           check=TRUE, run = TRUE, hashfilename=FALSE)

m0_model <- readModels("cfa_m1SB_uli.out")
output <- data.frame(m0_model$output) %>% slice(210:229) %>% rename(svalues=1) %>% summarise(svalues=paste(svalues,collapse="\n"))

#run new model with nested model parameters
cfa_m3SB <- mplusObject(
  
  TITLE = "CFA - ULI MOD for SB",
  
  VARIABLE = "usevar = rum-cnsq;",
  
  ANALYSIS = "estimator = mlr;
              convergence=100000000;",
  
  MODEL = glue("{output$svalues}
                cnsq with lfe*0;"),
  
   OUTPUT = "tech5;",
  
  usevariables = colnames(cfa),
  rdata = cfa)

m10_model_fit <- mplusModeler(cfa_m3SB,
                           dataout=("cfa_dataSB.dat"),
                           modelout=("m10_model.inp"),
                           check=TRUE, run = TRUE, hashfilename=FALSE)


#compute chi-sq test
M0 <- readModels("cfa_m1SB_uli.out")
M1 <- readModels("cfa_m3_uli.out")
M10 <- readModels("m10_model.out")

F0 <- M0[["summaries"]][["ChiSqM_Value"]]
F1 <- M1[["summaries"]][["ChiSqM_Value"]]

c0 <- M0[["summaries"]][["ChiSqM_ScalingCorrection"]]
c1 <- M1[["summaries"]][["ChiSqM_ScalingCorrection"]]
c10 <- M10[["summaries"]][["ChiSqM_ScalingCorrection"]]

d0 <- M0[["summaries"]][["ChiSqM_DF"]]
d1 <- M1[["summaries"]][["ChiSqM_DF"]]

f <- ((F0*c0 - F1*c1)/(d0-d1))/(c0*d0-c10*d1)

df_delta <- d0-d1

p_val <- pchisq(f,df_delta,lower.tail=FALSE)
p_val
```


### Residuals
```{r}
norm_res <- as.data.frame(out_mod2$residuals$covarianceResid.norm) %>% rownames_to_column("Variable") %>% replace(is.na(.), "")

norm_res %>% gt() %>% tab_header(
  title = md("**Table 6**"), subtitle=md("*Normalized Residuals*")) %>% cols_align(align="center") %>% 
  tab_options(column_labels.font.weight = "bold")


norm_res2 <- as.data.frame(out_mod$residuals$covarianceResid.norm) %>% rownames_to_column("Variable") %>% replace(is.na(.), "")

norm_res2 %>% gt() %>% tab_header(
  title = md("**Table 6**"), subtitle=md("*Normalized Residuals*")) %>% cols_align(align="center") %>% 
  tab_options(column_labels.font.weight = "bold")

norm_res3 <- as.data.frame(out_uli$residuals$covarianceResid.norm) %>% rownames_to_column("Variable") %>% replace(is.na(.), "")

norm_res3 %>% gt() %>% tab_header(
  title = md("**Table 6**"), subtitle=md("*Normalized Residuals*")) %>% cols_align(align="center") %>% 
  tab_options(column_labels.font.weight = "bold")
```

### Communality estimates (prop of variance in variable explained by factor)
```{r}
com_est <- as.data.frame(out_mod2$parameters$stdyx.standardized) %>% filter(grepl("Residual.Variances", paramHeader)) %>% 
  select(param, est) %>% mutate(est=1-est)

com_est %>% gt() %>% tab_header(title=md("**Table 8**"),
                                subtitle=md("*Communality Estimates*")) %>% cols_label(param="Variable",
                              est="Communality") %>% 
  tab_options(column_labels.font.weight="bold")

com_est2 <- as.data.frame(out_uli$parameters$stdyx.standardized) %>% filter(grepl("Residual.Variances", paramHeader)) %>% 
  select(param, est) %>% mutate(est=1-est)

com_est2 %>% gt() %>% tab_header(title=md("**Table 8**"),
                                subtitle=md("*Communality Estimates*")) %>% cols_label(param="Variable",
                              est="Communality") %>% 
  tab_options(column_labels.font.weight="bold")
```

### HIgher order factor
```{r}
cfa_h1 <- mplusObject(
  
  TITLE = "CFA - H1",
  
  VARIABLE = "usevar = rum-cnsq;",
  
  ANALYSIS = "estimator = mlr;",
  
  MODEL = "factor1 by rum strs ef dif dstr adj lie avd;
  factor2 by dstn intr aut;
  factor3 by obl exp glt;
  factor4 by uncf rep lfe cnsq;
  factor5 by factor1* factor3;
  factor6 by factor2* factor4;
 
  factor1@1;
  factor2@1;
  factor3@1;
  factor4@1;
  factor5@1;
  factor6@1;",
  
  PLOT = "type = plot3;",
  
  OUTPUT = "sampstat standardized residual modindices (3.84);",
  
  usevariables = colnames(cfa),
  rdata = cfa)

cfah1_fit <- mplusModeler(cfa_h1,
                           dataout=("cfa_data.dat"),
                           modelout=("cfa_h1_uli.inp"),
                           check=TRUE, run = TRUE, hashfilename=FALSE)


out_h1 <- readModels("cfa_h1_uli.out")
semPaths(out_h1, intercepts=FALSE)

model_fit_h1 <- LatexSummaryTable(out_h1,
       keepCols=c("Title", "Parameters","LL", "ChiSqM_Value", "ChiSqM_DF", "ChiSqM_PValue", "RMSEA_Estimate", "RMSEA_90CI_LB", "RMSEA_90CI_UB", "CFI", "TLI", "SRMR")) %>% 
  mutate_at(vars(contains("RMSEA")), ~format(., nsmall=3)) %>% unite(CI, RMSEA_90CI_LB:RMSEA_90CI_UB, sep = ", ", remove=TRUE) %>% mutate(CI = paste0("(", CI, ")")) %>% unite(RMSEA, RMSEA_Estimate:CI, sep=" ", remove=TRUE)

model_fit_h1 %>% gt() %>% tab_header(
  title = md("**Table 1**"),
  subtitle=md("*SUmmary of Model Fit Indices*")) %>% 
  cols_label(Title = "Model", Parameters = md("Par"),
             LL = md("*LL*"), ChiSqM_Value = md("CHi^2"),
             ChiSqM_PValue = md("*p-value*"), 
             ChiSqM_DF = md("*df*"),
             RMSEA = "RMSEA (90% CI)") %>% 
  tab_options(column_labels.font.weight = "bold") %>% 
  fmt(c(6), fns=function(x) ifelse(x<.001, "<.001", scales::number(x,accuracy=.01)))

varh1 <- out_h1$parameters$stdyx.standardized %>% 
  filter(grepl(".BY", paramHeader)) %>% 
  select(param, est, se)

fach1 <- out_h1$parameters$stdyx.standardized %>% 
  filter(grepl(".WITH",paramHeader)) %>% 
  select(param, paramHeader, est, se) %>% 
  unite(param, paramHeader, param) %>% 
  mutate(param=str_replace(param, "FACTOR6.WITH_FACTOR5", "F5 with F6"))

loadings_stdyxh1 <- full_join(varh1,fach1)

loadings_stdyxh1 %>% gt() %>% tab_header(
  title = md("**Table 2**"),
  subtitle=md("*Standardized Factor Loadings and Factor Correlation Estimates*")) %>% 
  cols_label(param = "Item", est = "Loading", se = "SE") %>% 
  tab_row_group(label = "Factor 1", rows = 1:8) %>% 
  tab_row_group(label = "Factor 2", rows = 9:11) %>% 
  tab_row_group(label = "Factor 3", rows = 12:14) %>% 
  tab_row_group(label = "Factor 4", rows = 15:18) %>% 
  tab_row_group(label = "Daily Demands", rows = 19:20) %>% 
  tab_row_group(label = "Anticipated Harm", rows=21:22) %>% 
  tab_row_group(label = "Factor Correlation", rows=23) %>% 
  row_group_order(groups = c("Factor 1", "Factor 2", "Factor 3", "Factor 4", "Daily Demands", "Anticipated Harm", "Factor Correlation")) %>% 
  tab_options(column_labels.font.weight = "bold")


cfa_h2 <- mplusObject(
  
  TITLE = "CFA - H2",
  
  VARIABLE = "usevar = rum-cnsq;",
  
  ANALYSIS = "estimator = mlr;",
  
  MODEL = "factor1 by rum strs ef dif dstr adj lie avd;
  factor2 by dstn intr aut;
  factor3 by obl exp glt;
  factor4 by uncf rep lfe cnsq;
  factor5 by factor1* factor4;
  factor6 by factor2* factor3;
  factor1 with factor4;
  factor2 with factor3;
 
  factor1@1;
  factor2@1;
  factor3@1;
  factor4@1;
  factor5@1;
  factor6@1;",
  
  PLOT = "type = plot3;",
  
  OUTPUT = "sampstat standardized residual modindices (3.84);",
  
  usevariables = colnames(cfa),
  rdata = cfa)

cfah2_fit <- mplusModeler(cfa_h2,
                           dataout=("cfa_data.dat"),
                           modelout=("cfa_h2_uli.inp"),
                           check=TRUE, run = TRUE, hashfilename=FALSE)


out_h2 <- readModels("cfa_h2_uli.out")
semPaths(out_h2, intercepts=FALSE)

model_fit_h2 <- LatexSummaryTable(out_h2,
       keepCols=c("Title", "Parameters","LL", "ChiSqM_Value", "ChiSqM_DF", "ChiSqM_PValue", "RMSEA_Estimate", "RMSEA_90CI_LB", "RMSEA_90CI_UB", "CFI", "TLI", "SRMR")) %>% 
  mutate_at(vars(contains("RMSEA")), ~format(., nsmall=3)) %>% unite(CI, RMSEA_90CI_LB:RMSEA_90CI_UB, sep = ", ", remove=TRUE) %>% mutate(CI = paste0("(", CI, ")")) %>% unite(RMSEA, RMSEA_Estimate:CI, sep=" ", remove=TRUE)

model_fit_h2 %>% gt() %>% tab_header(
  title = md("**Table 1**"),
  subtitle=md("*SUmmary of Model Fit Indices*")) %>% 
  cols_label(Title = "Model", Parameters = md("Par"),
             LL = md("*LL*"), ChiSqM_Value = md("CHi^2"),
             ChiSqM_PValue = md("*p-value*"), 
             ChiSqM_DF = md("*df*"),
             RMSEA = "RMSEA (90% CI)") %>% 
  tab_options(column_labels.font.weight = "bold") %>% 
  fmt(c(6), fns=function(x) ifelse(x<.001, "<.001", scales::number(x,accuracy=.01)))

varh2 <- out_h2$parameters$stdyx.standardized %>% 
  filter(grepl(".BY", paramHeader)) %>% 
  select(param, est, se)

fach2 <- out_h2$parameters$stdyx.standardized %>% 
  filter(grepl(".WITH",paramHeader)) %>% 
  select(param, paramHeader, est, se) %>% 
  unite(param, paramHeader, param) %>% 
  mutate(param=str_replace(param, "FACTOR6.WITH_FACTOR5", "F5 with F6"))

loadings_stdyxh2 <- full_join(varh2,fach2)

loadings_stdyxh2 %>% gt() %>% tab_header(
  title = md("**Table 2**"),
  subtitle=md("*Standardized Factor Loadings and Factor Correlation Estimates*")) %>% 
  cols_label(param = "Item", est = "Loading", se = "SE") %>% 
  tab_row_group(label = "Factor 1", rows = 1:8) %>% 
  tab_row_group(label = "Factor 2", rows = 9:11) %>% 
  tab_row_group(label = "Factor 3", rows = 12:14) %>% 
  tab_row_group(label = "Factor 4", rows = 15:18) %>% 
  tab_row_group(label = "Daily Demands", rows = 19:20) %>% 
  tab_row_group(label = "Anticipated Harm", rows=21:22) %>% 
  tab_row_group(label = "Factor Correlation", rows=23) %>% 
  row_group_order(groups = c("Factor 1", "Factor 2", "Factor 3", "Factor 4", "Daily Demands", "Anticipated Harm", "Factor Correlation")) %>% 
  tab_options(column_labels.font.weight = "bold")

mod_ind_h2 <- out_h2$mod_indices %>% arrange(desc(MI))

mod_ind_h2 %>% gt() %>% tab_header(title=md("**Table 3**"),
     subtitle = md("*Model Modification Indices*")) %>% 
  cols_label(modV1 = "Variable 1", operator = "Operator",
             modV2 = "Variable 2", EPC = "E.P.C.", 
             Std_EPC = "Std E.P.C.", 
             StdYX_EPC = "StdYX E.P.C.") %>% 
  tab_options(column_labels.font.weight = "bold")

```










